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ABSTRACT 

An analytic model is developed of baroclinic currents 
measured by several investigators under the Arctic Ocean ice 
cover. A time-dependent, Boussinesq, quasi-geostrophic 
vorticity equation is developed to describe the flow. 
Separable normal 'mode solutions are found which shov/ sta- 
tionary current patterns qualitatively consistent with 
observed data. An initial value problem is then formulated 
using density perturbations at a newly forming lead as a 
forcing function. 
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I. 



INTRODUCTION 



A. OBJECTIVF.S OF RESFARCH AND RFVIFW OF PRFVIOUS RFSFARCH 

The objective of this study was to develop an analytic 
model of specific small scale haroclinic current patterns 
which have been observed on a recurring basis beneath the 
Arctic ice cap. The patterns are termed small in that they 
are on a scale less than that of the general Arctic circula- 
tion. In order to develop the model a theoretical investi- 
gation was made of the baroclinic response of the Arctic 
Ocean to perturbations in its internal density structure at 
a newly forming lead. 

It is v-zorthwhile at this point to briefly discuss pre- 
vious research in the area of Arctic dynamic modeling, and 
to describe the baroclinic currents measured by several 
investigators, in an effort to indicate the relationship of 
this study to previous work. 

Investigations of Arctic Ocean dynamics can be classified 
by their scale as well as by the mode of motion described 
(either internal, surface, or a combined mode). The largest 
scale studied concerns processes which affect the entire 
ocean, such as long term mean ice drift and water circulation. 
In order to practically study motion on this scale, spacing 
between grid points cannot feasibly be less than tens of 
kilometers. This necessitates parameterizing some of the 
smaller scale processes occurring between grid points, rather 
than solving for them explicitly. Factors Vv’hich must be 
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considered on this scale are the shape and bathymetry of the 
ocean basin, Coriolis effects, wind and ice stress, etc. 

Several theoretical investigations of large scale Arctic 
Ocean dynamics which closely match observed patterns of water 
circulation and ice drift have been carried out in recent 
years. Among these, Campbell [1965] proposed a large scale 
analytic ice drift model and Galt [1973a] developed a working 
barotropic numerical model of overall surface water circula- 
tion. Both these models were concerned with surface motion 
or an equivalent barotropic circulation. A complete model 
of large scale internal Arctic circulation has not yet been 
published, but work along these lines is being carried out 
at NOAA's Geophysical Fluid Dynamics Laboratory, Princeton, 
New Jersey [Semtner, 1973]. 

The smallest scale of analysis concerns the short term 
local effects of phenomena such as refreezing leads and ice 
cracking. These effects result in dynamic processes such as 
forced convection due to the salt rejection and heat transfer 
during freezing. Grid size of concern here would be from 
less than a cm, a microscale, to a few meters. Coriolis 
effects would not be important at the smallest end of this 
scale but effects such as diffusion would have to be con- 
sidered. Several investigators have done work on the problem 
of forced convection associated with freezing of sea ice. 

Foster [1972] treated the problem analytically after 
having published several earlier papers on the problem of 
forced convection in the ocean. The initial state in his 



9 



study was a homogeneous under-ice layer. He then assumed 
density to be a function of salinity only and solved the 
equations of motion and salt diffusion directly associated 
with a refreezing lead. Using a two-dimensional model Foster 
carried his solution to the point of finding initial vertical 
velocities and length and depth scales for convection cells 
during the initial stages of convection. These were about 
1.4 cm/sec for vertical velocities, 80 cm for horizontal 
convection cell spacing, and several cms for depth of con- 
vective penetration. He did not analytically consider 
further stages of convection nor did he continue to a larger 
scale in which Coriolis forces would become important. He 
did postulate a higher order of convection cells with depth 
penetration of up to 6 meters and vertical velocities of 
about 1 cm/sec. 

Piascek [1970] considered the general problem of surface 
cooling effects on the ocean. He developed a numerical 
model for convection cells on both a small scale (order of 
cms) v;hich neglected Coriolis effects and a larger scale 
(order of kms) which considered Coriolis. His model con- 
figuration was a rectangular box directly below the surface 
with an initially homogeneous structure. He considered 
uniform surface cooling as a boundary condition and then 
developed a vorticity equation in terms of a vector poten- 
tial. This was numerically integrated over a 2-dimensional 
grid to give steady state temperature and vorticity patterns 
v;ithin his "box." With some modification his model could be 
adapted to describe Arctic under-ice convection. 
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Coachman [1966] considered a mass and energy balance on 
the layer of water iitimediately below the ice cover. He 
assumed a downv/ard flow direction under the refreezing leads 
due to sinking of highly saline v;ater and a compensating 
upward flow spread out over the area between leads , as shown 
in Figure 1. His study used mean values of freezing rate, 
temperatures, lead size and distance • between leads to deter- 
mine mean vertical velocity values. His scale of analysis 
was necessarily larger than that of Foster in order to arrive 
at mean velocities . His downward vertical velocity values 
are on the order of .01 to .08 cm./sec. The depth of his 
surface circulation layer was 100 meters. If the Arctic 
Ocean is viewed from a large scale standpoint, the net 
result of the investigations by Coachman, Foster, and Piascek 
would be to show that there is a small downward vertical 
velocity associated with each refreezing lead due to forced 
convection and that there is a local increase in the density 
of the v;ater directly under a refreezing lead. 

A related small scale problem which has been investigated 
is that of the ice formation rate and rates of salt rejection 
and heat transfer at a refreezing Arctic lead. Schaus [1971] 
developed a numerical model of a refreezing lead which treats 
these problems and which could be used in coupling the air- 
sea interaction problem with ice dynamics and v^?ith the under- 
ice convection problem discussed here. 

A vertical velocity similar to that caused by the convec- 
tive process may also arise from the initial cracking and 
ice motion during lead formation. This specific problem 
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Figure 1. 



Schematic of Convection Nodel Proposed 
by Coachman [19f6]. 
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has not been previously studied. In this case the theory is 
as follows, assuming ice shears upon cracking over a period 
of a few days as shown in Figure 2. There is a stress 
imparted to the underlying water at the ice-water boundary, 
which in turn causes a convergence of the Fkman layer. This 
in turn results in an impulsive vertical velocity under the 
newly forming lead. The net result is that, for an entirely 
different reason, a vertical velocity similar to that caused 
by convection is associated with shear across newly forming 
leads. Both these lead-associated phenomena would be classed 
as small scale effects. 

By interacting with the existing density and flow 
patterns these velocity and density perturbations may have 
some, as yet undetermined, effect on the large scale internal 
flow dynamics of the Arctic. A proposed theory for this is 
discussed and developed in part B of this section. Any 
resultant flow pattern due to such an interaction might 
then be expected on a scale intermediate in size between 
the smallest and largest scales discussed; on the order of 
meters in depth and meters to kilometers in length. A theo- 
retical investigation of baroclinic internal flow in the 
Arctic on this scale has not previously appeared in the 
literature. It is important to note here that this is the 
scale of a number of internal current patterns or events 
measured in the Arctic. 

The observed currents appear to be baroclinic in nature 
although a satisfactory theory for their origin has not yet 
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Figure 2. Proposed Shear Mechanism 
for Lead Forination. 
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been developed. An example of the vertical distribution of 
horizontal velocities measured by Galt [1967] under Ice 
Island T-3 is shown in Figure 3. The corresponding density 
distribution is also shown. Similar results found by 
Coachman and Newton [1973] are shovm in Figure 4. It can be 
seen that a peak in the horizontal velocity occurs in the 
pycnocline, vjith a maximum at about 150 meters, indicating 
a baroclinic origin. Figure 5 shows the horizontal distri- 
bution of velocities found by Coachman and Newton at 150 
meters. It is a plot of velocity vectors measured from 
one position on the ice while drifting over a given area. 

As can be seen, there is a periodic change in the velocity 
direction with an apparent wavelength of about 30 km. The 
reason for classing this periodicity as spatial rather than 
temporal is that the observed currents were found to be 
essentially in geostrophic balance. For this type feature, 
phase propagation speeds will be small compared to advective 
speeds. Figure 6 shows a vertical cross-section of veloci- 
ties as reported by Coachman and Newton [1973]. 

A need thus exists for an investigation of the Arctic 
Ocean internal flow dynamics on the intermediate scale 
described. The immediate purposes of this study are to offer 
an explanation of the observed baroclinic currents. A 
further use of such a study is to couple the small scale 
phenomena to the overall circulation, which in turn may aid 
in the development of a three dimensional Arctic circulation 
model . 
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Figure 3. Vertical Distribution of Horizontal 
Velocities Measured by Galt [1967]. 
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Figure 4. Vertical Distribution of Horizontal 
Velocities Measured by Coachman and 
Nev/ton [197 3] . 
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Figure 5. Horizontal Distribution of Horizontal Velocities 
Measured by Coachinan and Mewton [1973] . 
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Figure 6. 



Vertical Section of Velocities Measured 
by Coachman and Nev/ton [1973]. 
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B. THEORY DEVELOPED IN THIS THESIS 



A theory was postulated to describe the dynamic effects 
on the underlying ocean of a nev/ly forming Arctic lead. The 
theory developed was begun with the premise that a small 
vertical velocity perturbation or density increase is created 
during the formation of each lead, caused by either refreez- 
ing or cracking, as described earlier. The scale of analysis 
used v/as tens of kms in the horizontal and the full depth 
of water (0-3000 m) in the vertical. 

The water directly below the Arctic ice cover can typi- 
cally be expected to consist of a shallow homogeneous layer 
(20-50 m) , under which a stably stratified density distribu- 
tion exists. There is also some horizontal slope to the 
isopycnals in order to support a mean flow. The horizontal 
gradient will be several orders of magnitude less than that 
of the vertical density gradient in the pycnocline. Pertur- 
bations such as the vertical velocity and density changes 
occurring under a lead, which is a very small area on the 
assumed length scale, should distort the isopycnals over a 
relatively short period of time and short distance. This 
would act as an impulsive event on an overall time scale 
measured in days. As the isopycnals attempt to return to 
their stable condition, a transient waveform of some sort 
would then appear and propagate out from the lead. 

The first problem would then be to find the form of the 
waves which can exist in a typical Arctic density structure 
and to then determine if these are consistent with observed 
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current patterns. Once this was done, a second problem would 
be an investigation to determine whether the perturbations 
at a lead could actually force these oscillations. An inves- 
tigation of these two problems forms the basis of this thesis. 

An indication of the possible validity of the work at any 
point is the degree of correlation of theoretically developed 
waves or currents v;ith real data. A comparison was therefore 
made with the baroclinic currents measured by Galt [1967] and 
Coachman and Newton [1973] once a theoretical current pattern 
was developed. 

The first step in substantiating the theory was to inves- 
tigate the type of waves which could exist in the given 
density structure and to compare these with known data. This 
corresponded to treating the Arctic Ocean as a semi-enclosed 
system, subject to certain governing equations of motion and 
boundary conditions, and solving for its preferred modes of 
oscillation. The actual solution to this problem is the 
subject of Section II of this thesis. The general assumptions 
made in formulating the theory and the rationale behind them 
will be discussed here. 

The equations used follow the development of Galt [1973b] 
which is based on that of Kuo [1973]. The basic difference 
from Kuo's development is in the formulation of the equation 
of state. The system was assumed to be inviscid, hydrostatic, 
and quasi-geostrophic . The equation of state assumes density 
proportional to salinity only. This is a reasonable approx- 
imation in the region being investigated as can be seen in 
Figure 7, which is a plot of salinity, temperature, and 
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Figure 7. Observations of Salinity, Temperature and 
Density Taken Under Ice Island T-3. 
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density profiles taken under T-3 along with current measure- 
ments [Galt, 1967]. 

The Boussinesq approximation was also made, V7hich states 
that density variations are significant only when the density 
is multiplied by gravity. An additional assumption necessary 
to this development is tliat the Rossby number is small. This 
allowed the expansion of the dependent variables as power 
series in terms of the Rossby number. Equations consistent 
in powers of Rossby num}:)er could then be derived. 

From the above assumptions and the basic equations of 
motion a quasi-geostrophic vorticity equation was developed 
in terms of the perturbation pressure (pressure could be 
written as the sum of a mean and a perturbation quantity) . 
This equation was solved by assuming a separable solution 
consisting of the product of a horizontal plane wave and an 
unknown vertical structure function. This vertical structure 
function was then solved for numerically s\ibject to appro- 
priate boundary conditions. Results of the numerical solu- 
tion are compared with the measured currents in Section III 
of this thesis. 

Investigation of the second problem described (whether 
the theoretical oscillations could actually be driven by 
perturbations found at a lead) was accomplished by using 
pressure conditions under a newly forming lead as an initial 
condition for the quasi-geostrophic vorticity equation 
developed in Section II. A description of this initial value 
problem forms Section IV of this thesis. 
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II. 



DF.VELOPMFNT OF THE WAVE SOLUTION 



A. ANALYTIC DEVELOPMENT OF THE VORTICITY EQUATION 



The initial formulation of the problem was carried out 
using a rectangular coordinate system, with positive values 
of X to the East, y to the North, and z upward. The govern 
ing equations of motion for the flow can be written in the 
following Boussinesq, inviscid form: 



u. + uu + vu + wu 
t X y z 



V. + uv + w + wv 
t X y z 



+ UW + W7 + ww 
t X y z 




( 1 ) 

( 2 ) 

(3) 



where u,v, and w are the x,y, and z velocity components, 
respectively, f is the vertical component of the Coriolis 
parameter, g is the acceleration of gravity, p "is the pres- 
sure, p is the density and is the horizontal space and 
time average of the density, which is a function of z only. 
In addition to the equations of motion, the continuity equa 
tion can be written as : 



p. + up + vp + wp + p (u +v + w ) = 0 (4) 

^t *^x ^y ^z ^s X y z 

and the equation for conservation of salt, assuming a con- 
servative lav; with no mixing, is: 

s. + us + vs + ws =0 (5) 

t X y z 
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Equations (1) through (5) can then be scaled and non- 
dimensionalized using essentially the same system as was 
used by Kuo [1973]. The non-dimensional quantities (primed) 
are defined as follov7s : 

(x,y) = L(x',y'); (u,v) = U(u',v'); (6) 




P = Pg + f^ULp^cf)' ; p = pg (1 + ep' ) ; 

s = s ( 1 + es ' ) ; 
s 

where L is a characteristic horizontal scale, D is a charac- 
teristic depth scale, U is a characteristic horizontal 
velocity and p^ and s^ are the static pressure and density 
fields respectively. 

In addition the following non-dimensional coefficients 
are defined; 

R = g — , Rossby number; (7) 

o'^ 

fo^L^ 

y = ^ , a ratio of the scale length, L, to 

^ the Rossby deformation lengtli; 

e = yR , a measure of the perturbations in 
the thermodynamic quantities; 



5 

P 



JL s 

— — , a measure of the stability of the 

*^s ^ static density field; 



<S 

s 



s 

s 



, a measure of the stability of the 
static salinity field. 



Substitution from (6) into (1) and (2) then gives: 
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( 8 ) 



R(u^, 


+ u * 


u' , + 


V ' u ' , + Rv/ 

y 


■u;.) 


I 

< 






U * 


v' , + 

X ' 


V ' V ' , + Rw 

y 


-v;.) 


+ f'u' 


= -4>' . 
y 


and substitution 


from 


(7) into 


(4) gives: 








, + u ' 


w' , + V ' w ' 

X ' y 


, + Rw ' 


'w;.)= 


-(j)' , - 6 



(9) 



providing that the hydrostatic relationship for the steady 
state density and pressure fields is used, as given by: 



'Psh, “ 



Then, substitution of (6) into (4) gives: 



e ( Pt I + u ' p', + V ' p' , + Rw ' p' , ) + — (w'p ) ,+u', + v', = 0 (11) 

^t' y' z' Pg "^s z X y' 



Finally, substitution from (6) into (5) gives: 



R(s' +u's', + v's', + Rv;'s', + 
t X ' y ' z ' 



■fl ' 



( 12 ) 



Equations (8) through (12) can then be simplified by consi- 
dering limitations which occur in tlie real ocean: 



( — 1 

V 

V 


used 


in eq . 


(10) 


(13) 


e = pR « R 


used 


in eq. 


(11) 


(14) 



and by assuming a simplified equation of state relating 
density to salinity: 

p = s 

or, in a non-dimensional form for the perturbation quantities: 

p' = s' ; 6 = 6 = 0 (R) . (15) 

P s 
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Using eq. (15), the non-dimensional equations become 
(omitting primes from here on) : 



R(u^ 


+ uu + 

X 


vu 

y 


+ Rwu ) - fv 
z 


= ~(l> 

^X 


(16) 


R(v^ 


+ uv + 

X 


vv 

y 


+ Rwv ) + f u 
z 


II 

I 


(17) 


♦z 


6 (J) + s 
P 


= 0 




• 


(18) 


u + 

X 


+ 

>1 

> 


(wp^ 


N 

II 

O 




(19) 














R(St 


+ us + 

X 


vs 

y 


+ Rws^ + |-jf| 


w) =0 


(20) 



At this point the assumption of a small Rossby number, R, 
will be used and dependent variables expanded in a power 
series form; i.e., u, v, and s will take the form: 

a = E R^a (21) 

n=0 

where the a represents the general dependent variable and w 
will be written as : 

2 

w = w^ + Rw 2 + R w^ + ... etc. 

This is consistent with the scaling given by equation (6) . 
In addition, the Coriolis parameter will be scaled on a 
beta-plane : 

f^ = 1 + (22) 

where 

By £ 0(R) (23) 

Using these approximations plus eq . (15) and the series 
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expansion gives the following equations to zero order in 
Rossby number: 



-V 

o 




u 

o 



= -(j) 



o 



y 






o 



z 



-s 



o 



U + V 

o o 



0 



(24) 

(25) 

(26) 
(27) 



X y 



These show the flow to be geostrophic, hydrostatic, and hori- 
zontally non-divergent to a zeroth order approximation. The 
equations to first order in Rossby number are: 



u +UU +VU -Vt~ Byv = -(f), 

o. o o o o 1 o ^1 

t X y 



X 



V +UV +VV +u,+ 3yu = -d), 
o. o o o o 1 o ^ ly 

t X y 



= -6 (b - s, 
^1 p^o 1 

z 



^1 '"l ® 

X y s ^ 



S +US + VS + 

o. o o o o 

t X y 



I vv 

h / 1 



w, = 0 



(28) 

(29) 

(30) 

(31) 

(32) 



To eliminate the first order pressure terms from (28) 
and (29) the equations are cross-differentiated and sub- 
tracted, resulting in: 



(v -u ) + u (v -u ) + V (v -u ) + (v -u ) (u +v ) 

o,, o t o o o o o o o o o o 

xy^ xy xy .xyxy 



+ (u, +v, ) + 3y (u +v ) + 3v =0 

1.. 1_. o o o 



(33) 



X y 



X y 



Then, substitution from (27) and (31) into (33) results in: 
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( 34 ) 



(v -u ) + u (v -u ) -I- V (v -u ) 

O O t o o o o o o 

X y ^ X y X y 



- — (w, p ) + 3 v = 0 



p ' I's 
s 



o 



Equation (32) is rewritten in the follov;ing form; 



(WjP^) + 






(s +US +VS )=0 
o. o o o o 

t X y 



(35) 



Substitution of (26) into (35) (noting that -g 7 ^ f (x,y,t)) 

'^P 

gives : 



(Wj^p^) = 



0 + u 

0 o , o 

’ p z' t 



S ^o 



+ V 



pp. 



<f>o I (36) 



X ' P Z' y 



and differentiating this with respect to z gives: 






yp 

p ^z 



S A 

6 — 'f'o 



+ u 



zt 



^Ps 

S o o 

P Z zx 



pp. 



6 ^o / 

P z zy 



+ u 

o 

z 




o 



+ V 



o 





(37) 



which can then be substituted directly into (34) . In writing 
this it is helpful to define all of the velocities in terms 
of the pressure, assuming a geostrophic relationship (this 
assumption filters out gravity waves) . 



u 

o 




V = (j) 

o ^o 

X 



(v -u ) = (b +6 
o o XX ^yy 

X y 



(38) 



From these it can be seen that the last two terms in (37) 
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identically cancel. Then, substitution of (37) and (38) 
into (34) and dropping all zero subscripts gives: 



q^ + (f)q - (i)q = 0 
^t x^y ^y^x 



( 39 ) 



where : 




z 



(40) 



Equation (39) is the quasi-geos trophic vorticity equa- 
tion and (40) is the potential vorticity. The first two 
terms on the right hand side of (40) constitute the relative 
vorticity. The third terra is the contribution due to the 
variation in the Coriolis parameter, and the last term repre- 
sents the vorticity contributed by vertical stretching of the 
water column. Equation (39) shows that time-dependent motion 
will result when the motion is such as to create a change in 
potential vorticity. For the initial problem studied here, 
the Coriolis parameter, 3, was assumed to be non-zero and 
the resultant motion investigated. It should be noted that 
as the North pole is approached, 3 approaches zero, but that 
motion could still be driven by bathymetric variations. 

Assuming a separable v;ave solution of the form: 



leads to the following Sturm-Liouville problem for the verti- 
cal dependence: 






(41) 




( 42 ) 
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Substituting from eq. (7) for y and 6 , using the fact 

p T 9p ■ 

the Brunt-Vaisala frequency, N, is given by - 

P o Z 

converting back to a dimensional form results in : 



that 

and 




z 

z 

P 




k2 + i2 + Bk 


s 2 

L s 




a ^ 



Z 



0 



(43) 



B. EQUATION USING ARCTIC DENSITY PROFILE 
1 . Basic State Description 

In this initial formulation zero mean flow was 
assumed, thus the initial density distribution was assumed 
to be a function of depth only. The Arctic under-ice density 
profile was approximated by a homogeneous upper layer from 
the surface to z = -h , belov? which density increased follow- 
ing an exponential curve. In this lower layer density 
increases such that it assymptotically approaches a constant, 
, at deeper depths. This distribution is similar to that 
proposed by Reid [1948] , in his equatorial model, and is 
given analytically by: 



= P 



h 



(0 > z > -h) 



(44) 



Po = P,q - 



'Od-fh* 



1/c (z+h) 



(z < -h) 



which is illustrated in Figure 8 for various values of c and 
h (referred to hereafter as layer depth (h) and e-fold depth 
(c) ) . Actual density values taken during the collection of 
velocity data are also shown on Figure 8 for comparison. As 
can be seen, the best data fit is v;ith an e-fold depth of 
botv/een 100 and 120 meters and a layer depth of 30 meters. 
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Figure 8. Depth Distribution of Analytic Density 
Function Shov.’ing Relation to Data. 
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2 . Use of Actual Density Distribution in Fquations 

The density distribution given by eq. (44) was sub- 
stituted into eq. (43) to produce a vertical normal mode 
equation which would apply to this specific density structure. 
To shov; this, eq. (43) is first expanded in terms of and 
its derivatives: 



f 

o 



2 



g 




+ 2Z 



z 



p p 
s s 



zz 



■(Ps ^ 

z 



ZJ 



+ 




.+ 1 ‘ 




(45) 



The derivatives of are found from eq. (44) as: 

“^Pd“Ph^ l/c(z+h) 
p = e 

^z ^ 

“(Pd“Ph) l/c(z+h) 



(46) 



which can then be substituted into eq. (44) to yield, after 
some rearrangement of terms: 



Z 



zz 



^^Pd Ph^ l/c(z+h) ^ 1 

pc CJ 

V S 



f c 
o 



r <pd-0h> 



+ d + » 

a 



_^l/o(z+h)^ 



(•17) 



This can be simplified to: 



zz 



Ael/c(=i+h) + B 



Z + Xgl/c (z+h) 2=0 (48) 



where : 



A 



^ ^Pd~Ph^ 



P c 
s 



(49) 



33 



B 



1 

c 



(49) 



X 




'fd-Ph’ 




The solutions to eq. (48) will then be a set of eigen- 
functions describing the vertical distribution of the pressure 
function, (f) . It is noted that this equation only applies to 
the region below the layer depth. In the upper homogeneous 
region horizontal velocities will be constant since there is 
no stratification. Vertical boundary conditions on the 
system are that vj, the vertical velocity, is equal to zero at 
both the surface (ice) and bottom boundaries. In addition, 
velocities are assumed to be continuous across the internal 
boundary at z = -h. If these boundary conditions are met 
and a solution to eq. (48) exists there will be solutions 
for different discrete values of X, the eigenvalues for this 
equation. 

A value of either k or o may then be assumed and used 
with the eigenvalues to solve for the corresponding values 
of a or k for each mode. In order to actually solve for both 
k and a some initial distribution for (j) must be known and 
matched. This is the initial value problem, based on the 
perturbation at a lead, which is the next phase of the study 
and is discussed in Section IV. 

The problem being discussed here was formulated with 
the assumption of a long narrow lead aligned in a North-South 
direction. Thus x^ « y^ , where x^^ and y^ are the Fast-Nest 
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and North-South dimensions of the lead, respectively. Based 
on this it was then assumed that no changes occur in the y 
direction (9/3y = 0) which in turn means that 1 is assumed 
to be zero in eqs . (41) and (43) and that u =-3({)/3y = 0. 

The equations thus become two-dimensional. A horizontal 
boundary condition at time zero of v = 0 at the origin, 
which is placed in the center of the lead, is then imposed. 
The horizontal velocity can be written as: 



M. 

3 



ikZ (z) e 



i 



(kx-CTt) 



or : 



V = kZ ( z) [C^cos (kx~at ) + C 2 sin (kx-at ) ] 



For v(0,0) = 0, must equal zero or: 

V = kZ (z) sin (kx-at) (50) 

This means that (j> can be written as : 

(j) = Z ( z ) cos (kx-at ) 

From eq. (36) w can be written as: 

^Ps 

w = -? Z sin (kx-at) 

0 z 
P 

which becomes : 
f 

w = — ^ Z sin (kx-at) (52) 

N^ 2 



when the terms are dimensionalized and proper substitutions 
made for y and 6p. Assuming a value for k is known, an 
expression for a can then be derived from eq. (49) : 
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a 



( 53 ) 



-Bk 



f CD X 

o s 

g(P<j-Ph> 



+ 



This is similar to the equation for the frequency of a baro- 
tropic Rossby wave, = 3/k , with the additional term in 
the denominator being due to the stratification. 
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III. RESULTS OF WAVE SOLUTION AND CO>'’PARISON 
WITH ACTUAL DATA 



A. NUMERICAL SOLUTION OF THE VERTICAL MODE EQUATION 

The vertical mode equation (48) was rewritten using a 
normalized depth variable, 0 ^ z' ^ -1.0, as: 

Z - (54) 

z z z 

2 

where X ' = D X and D = depth of water . 

Equation (54) was solved numerically on the Naval Post- 
graduate School IBM/360 computer. The solution was found 
using a Milne predictor-corrector method, according to an 
algorithm described by Nielsen [1964]. The depth range was 
divided into 2000 increments and then the Milne method used 
from the bottom, z' = -1.0, upward to the normalized layer 
depth, z' = - h/D. The Milne method section of the numerical 
solution v/as then tested by using it to solve a zero order 
Bessel equation which had a known analytic solution. Results 
of this test were satisfactorily checked with a table of 
Bessel functions. 

In the relatively thin homogeneous region above the layer 
depth, the slope of the vertical velocity was. held at its 
layer depth value, as was the value of the vertical structure 
function, Z. The output of the numerical scheme was a set of 
2000 values of Z and Z^ from z' = -1.0 to z' = 0, for a given 
X'. The vertical distribution of horizontal velocity is 
proportional to Z as was shown by eq . (50). From eq . (52) it 
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can be seen that the vertical dependence of vertical velocity, 

2 

w, is proportional to Z /N . The numerical scheme could thus 

z 

be used to give the vertical dependence of both horizontal 
and vertical velocities . 

The bottom boundary condition, w = 0 at z ' = -1.0, was 

met by setting Z (-1.0) = 0. This was valid since w is 

z 

2 2 

proportional to Z /N and since N (-1’. 0), although quite 

z 

small, is not identically zero. With the bottom boundary 
condition thus fixed, the condition w = 0 at z ' =0 could 
only be met for certain values of X ' . These values are the 
eigenvalues corresponding to the normal modes of oscillation 
of the ocean with the given vertical density distribution. 

In order to find the correct eigenvalues an initial X' was 
assumed. An iterative process was then begun in which the 
numerical method described earlier was used to ari'ive at a 
value of w(0) for each new value of X'. The initial value 
chosen for X ' was then increased in steps until the value 
of. w(0) changed sign. At that point a Regula-Falsi method 
(eq. (55) below) was used to find the next X'. 



^n+1 



n-1 n 



w (0) - 



w (0 ) 
n 



X'w^ (0) 
n n-1 

w , (0) 
n-1 



(55) 



This was continued until the value of w(0) was sufficiently 
close to zero. The value of X ' at that point was then con- 
sidered an eigenvalue. The process was continued by increas- 
ing the value of X ' in steps until the next eigenvalue was 
found in the same manner. 
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B. RESULTS OF THEORETICAL DEVELOPMENT 



Figures 9 and 10 show the vertical distribution of 
horizontal and vertical velocities which result from the 
numerical solution of eg. (54) for a given density distribu- 
tion (h=30, c=120). The first four normal modes are shown 
in each figure. Velocity values plotted are not absolute 
values but are relative to the horizontal velocity at the 
bottom for each mode. The ratio of w(z) to V(D) was computed 
using a and k values corresponding to a wavelength of 30 km. 
As can be seen on these graphs, the vertical velocity would 
be insignificant for a horizontal velocity on the order of 
cm/sec or larger. This means the only measurable motion 
would be in the horizontal direction. 

Figure 11 shows the effect of a change in the density 
distribution layer depth on the horizontal velocity distribu- 
tion. The second mode only is shown since effects on other 
modes are similar. The values chosen, 20 and 50 meters, 
represent the normal seasonal variation in the thickness of 
the under-ice layer. As can be seen, very little change in 
the velocity distribution results. 

Figure 12 shows the effect on the second mode of a change 
in the curvature (e-fold depth, c) of the below layer density 
distribution. This shows that as the e-fold depth increases 
(decrease in curvature) , the velocity peak is less marked 
and occurs at a deeper depth. 

Frequencies v>;ere solved for by using eg. (53) v:ith an 
assumed horizontal wavelength of 30 km. Results were on the 
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Figure 9, Vertical Structure of Theoretical Horizontal 
Velocity Patterns for Modes 1 through 4. 



Yerfical Veloctfy v C3000^ 




Figure 10. Vertical Structure of Theoretical ^7ertical 
Velocity Patterns for Modes 1 through 4. 
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Figure 11. Fffect of Layer Depth Change on Vertical Structure. 



42 



Horiron«al Velocity v Cz)/ vC3000 ) } 



-20 -15 -10 -5 0 5 10 15 




Figure 12. Effect of Varying e-fold Depth in Density 
Distribution on Vertical Structure. 



~ 8 ~ 1 

order of 10 sec or less, which corresponds to periods of 
tens of years. These values indicate the current patterns 
are essentially stationary. As can be seen from eq. (53) 
either a decrease in wavelength (increase in k) or an increase 
in mode (higher X) will cause the frequency to decrease. 

Thus the limiting, highest possible frequency for a given 
density distribution and horizontal Wavenumber is for X = 0, 
a barotropic Rossby wave mode, and is given by: 




For a wavelength of 30 km this corresponds to a period of 
about 10 years. 

The horizontal distribution of velocities was arrived at 
by assuming a length of 15 km to the first zero of the 
sin(kx-at) function of eq. (50) and forming the product of 
this trigonometric function and the Z(z) function at a con- 
stant depth. This is shown in Figure 13. 

C. COMPARISON WITH OBSERVFD DATA 

A comparison of the theoretical normal modes plotted in 
Figure 9 with the observed data in Figure 2 shows that the 
second normal mode most nearly matches the measured vertical 
distribution of horizontal velocity. The only real discre- 
pancy between the two is the relative maximum in the theo- 
retical distribution at the surface. This appears to be an 
inherent characteristic of the theoretical solution and is 
seen in each of the normal modes. One possible explanation 
is that these theoretical maxima were caused by treatment of 
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Figure 13, Horizontal Distribvition of Theoretical 
Horizontal Velocities. 
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the under-ice surface as flat and frictionless in the develop- 
ment. The ice actually has a rough surface capable of dis- 
sipating some energy of relatively slow moving currents. 

A second explanation is that more than one mode exists 
with a progressively decreasing amount of energy in succes- 
sively higher modes. Thus if represents the horizontal 

velocity distribution for a given mode, then: 

OG 

V (z) = E a V (z) , 

, n n 
n=l 

where a > a , , , v/ould describe the actual velocity distri- 
n n+1 ^ 

button. Assuming further, as a simplistic first approach, 
that the only significant a^ are for n=l and n=2 (i.e., only 
the first two modes contain significant energy) the linear 
combination of modes 1 and 2 was plotted in Figure 14 . The 
a^^ and were sized to give zero velocity at the surface, 
meaning '^sls less than a^. Figure 14 can be seen to be 
closer to the observed data than the second mode alone. 
Obviously a more rigorous solution for the coefficients could 
be used to get an even closer fit to the data. This, however 
would not be physically meaningful. The purpose of this 
initial investigation was merely to show that the theoretical 
development would allow motion similar to what has been 
observed. This has been shown. The next phase of the inves- 
tigation, which will be discussed in Section IV of this 
thesis , is to study the initial value problem of how this 
specific ocean system responds to the initial perturbations 
in its density structure caused by a newly forming lead 
(i.e., which normal modes, if any, are excited). 
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Figure 14. Linear Sum of Modes 1 and 7 for 
h = 30, c = 120. 
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In order to compare real and theoretical distributions 
more closely, a vertical cross-section through the origin 
was drawn showing theoretical horizontal velocity values , 
Figure 15. This was done for the second normal mode for a 
density distribution having h=30 and c=120, and for the 
sim.ple horizontal distribution of Figure 13. Thus Figure 15 
could be compared with Figure 4, which is drawn for real 
current data having a similar density distribution. A 
similarity can be seen in the general flow direction and 
in the shape of the isolines. Apart from the relative 
smoothness of the theoretical plot, the only significant 
difference is in the upper homogeneous layer, v/hich has been 
discussed previously. Figure 16 is the same cross-section 
drawn for the summation of tlie first tvjo normal modes. As 
can be seen, the cell structure is closer to the observed 
data of Figure 6 than was the second mode alone. 
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Figure 15. Vertical Cross-section of Horizontal 
Velocities^ Mode 2. 
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Figure 16. Vertical Cross-section of Horizontal 
Velocities, Sum of Modes 1 and 2. 
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IV. 



DEVELOPMENT OF THE INITIAL VALUE PROBLEM 



A. GENERAL METPiOD 

The problem addressed here is whether the density field 
perturbations which occur at a newly forming Arctic lead can 
force the theoretically developed oscillations to occur and, 
if so, whether the forced oscillations are consistent with 
observed data. 

The solution to this problem lies in formulating expli- 
citly the initial effect of the lead-induced perturbations 
on the pressure field. This effect can then be used as an 
initial condition on the pressure function, cj) , in order to 
fully solve the Sturm-Liouville problem developed in Sections 
II and III. It v/as shown that the general solution to tliat 
problem v/as a set of eigenfunctions which gave the initial 
response of the Arctic Ocean to a baroclinic perturbation. 

The linear combination of these functions can be written as: 

oc oc 

<})(x,z,t) = E E A B Z (z)cos(k x-a t) (56) 

^ , n m n m mn 

n=l m=l 

so that at time zero: 

00 oc 

(j) (x,z) = E E A B Z (z)cos(k x) (57) 

^o , , n m n m ' 

n=l Pi=l 

It must be noted that the summation of solutions indicated by 
eq . (56) is not necessarily a solution itself. It would only 

constitute a solution if the initial equations v/ere actually 
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linear, which they are not. It appears reasonable to assume 

however that the non-linear contribution is small compared 

to the linear so that this summation approximates an actual 

solution. A method will nov/ be shown whereby the values of 

A ,B ,k and a can be found based on the initial specifica- 

tion of the ^ field. 

^o 



B. SPECIFICATION OF INITIAL CONDITIONS 

The perturbations at a lead during formation and refreez- 
ing were described as small downward velocities and an 
increase in density under the lead due to salt rejection and 
small dov7nward velocities due to shear processes. Consider- 
ing the increase in density below the lead as the dominant 
process, the initial density perturbation will result in a 
deviation from the basic state as shov/n in Figure 17; where 
Zp shown is the depth of penetration of the density pertur- 
bation. Assuming that this takes place only directly under 
the lead, then a vertical cross-section of initial density 
would appear as shown in Figure 18. Using the hydrostatic 
assumption, this will result in an initial overpressure con- 
dition (perturbation pressure p^ = as shovm in Figure 19 

for -X < X < X and zero elsev;here. This can then be used 
P - - P 

as an initial condition to which eq. (57) can be fitted in 
order to solve for the appropriate wavenumbers, frequencies, 

« 1 , will 



Z “’ll 

and constants. A further assumption that ' p 



be made, thus z^ - -h can be used in describing the initial 
distribution analytically. (})^(x,z) can thus be written: 
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Figure 17. Perturbation Density Due to Salt 
Rejection at Tirne Zero. 




X 



Figure 18. Vertical Cross-section of Perturbation 
Density at Time Zero. 
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Figure 19. Perturbation Pressure 
at Time Zero. 
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(f)^ (x,z) 



h 



X < X 



0 < z < z 

- p 



= 1 



1 

z > z 



(|x| > Xp) 



(58) 



C. SOLUTION OF THE PROBLEM 

The initial condition, <})^(x,z), is separable and can be 
written as (j)^(x,z) = f(x)g(z) v;here: 



f (x) =0 


(lx| 


> X ) 

i p 


(59) 


= 1 


(|x| 


i < X ) 

i - p' 





g(z) = -z/h (z > -h) 

=1 (z < -h) 



(60) 



Thus, coefficients can be solved for separately from: 



f (x) = E B cos (k x) 
, m m 

m=l 

oo 

g(z) = E A Z (z) 

, n n 
n=l 



(61) 



Considering B^^ first, these are the coefficients for a 
Fourier cosine series, which on the interval (0,L) are given 
by : 



B = 



B 

m 



X 

^ / f (x)dx + r / f (x)dx = 

^ 0 ^ X ^ 

x P 

^ ^ f (x) cos dx (m=l ,2,3...) 



(62) 



= (m=l,2,3...) 



(63) 
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where k = imr/L is a distance sufficiently far removed from 
m 

the lead so that (p can be assumed to have a zero value. 

The values of correspond to the zeroth mode, which is 
a barotropic mode that may exist but is not of interest in 
the present discussion. Figure 20 is a plot of the approx- 
imation to f (x) formed by evaluating eg. (61) for m = 1 
through m = 40. It can be seen to approximate f (x) fairly 
closely. Figure 21 is a plot of the horizontal distribution 
of V for these same modes . 

For each thus found, and each eigenvalue corresponding 

to the Z (z) vertical structure functions, there can be found 
n • 

a value of o 

mn 

The values of are given by: 

-D 

/ g(z)Z^(z)dz 

A = (n=l,2,3...) (04) 

/. [Z (z) ] ^dz 
0 

for a general series expansion, assuming the 7j^{z) are 
orthogonal functions with respect to a weighting function 
of one. The denominator can be normalized for all modes by 
finding a factor such that: 

2 

/ [C^Z^ (z) ] dz = 1 



or : 



/ (Z^(z))^dz 

0 

Thus A becomes : 




( 65 ) 
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Figure 20. Fourier Series Fxpansion of Initial 
Horizontal Pressure Distribution. 
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Figure 21. Fourier Series Fxpansion of Horizontal 
Velocity Distribution in ,x-direction . 
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( 66 ) 



-D 

A = / g(z)C Z (z)dz (n=l,2,3...) 
n Q n n 



for the normalized distributions. The computation of A^ can 
be further simplified by writing g(z) as the sum of an 
average and a variable part. Since D » h, the average part 
is essentially; 

g (z) - 1. 

^avg 



The only variable part is then: 



g (z) = 1 + z/h (z > -h) 
^var 



(67) 



Thus, eq. (64) can be written: 



-D 



-h 



A = / C Z (z)dz + / (l+z/h)C Z (z)dz 
n Q n n ^ n n 



( 68 ) 



But, the first integral in eq . (68) is equal to zero based on 

the original continuity relationship used in deriving Z (z) . 
From Figure 9 it can be seen that each Z (z) is constant for 
depths above -h. Thus eq . (68) finally becomes: 



A = - h/2 [a C ] (n=l,2,3...) 

n n n 



(69) 



These values were solved for using the IBM/360 computer and 
were then used to form the series summation given by eq. (61) 
for the first six modes. The results are shown in Figure 22. 
The approximation is not as close as the horizontal case due 
to using only six modes. This was a limit imposed by the 
numerical method. 

A complete set of values for A , B , k , and o is then 
available. Values for the distribution having c=120, h=30, 
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Figure 22. Series Fxpansion of Vertical 
Structure Function. 
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are given in Tables I, II, and III. Values of L=20 km and 

r =1 km were assximed. 

P 

The solution for ({) given by eq . (56) will then begin 
with a distribution as given by eq . (57), which is the ini- 
tial pressure distribution. Due to the extremely small 
values of a there will be no significant movement with time 

of the overall wave system. Phase velocities (a/k) are on 

_2 

the order of 10 cm/sec. The distribution of horizontal 
velocities however, is proportional to: 

00 oo 

v(x,z,t) = S E Vnt^)Bin(k x-a t) (70) 

n=l m=l 

The initial horizontal and vertical distributions of v(x,z,t) 
were shown in Figure 21 and Figure 22. 



TABLF I 

VALUFS OF (sec“^xlO“^) AND T^j^(yrs) FOR 
SELECTED VERTICAL AND HORIZONTAL MODES 



ml 5 10 



1 


9.83 


5.84 


3.11 


a 

mn 




20.0 


33.9 


63.4 


T 

mn 


n 3 


2.58 


4.38 


2.85 


a 

mn 




76.9 


45.2 


69.3 


T 

mn 


5 


1.11 


3.02 


2.49 


a 

mn 




178. 


65.6 


77.0 


T 

mn 
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TABLE II 



VALUES OF k AND B FOR HORIZONTAL 
m m 

MODES 1 THROUGH 10 



m 


km(cm ^xlO 


Bm 


1 


1.56 


. 0995 


2 


3.14 


.0983 


3 


4 . 70 


• .0963 


4 


6.28 


.0935 


5 


7.85 


.0900 


6 


9.40 


.0858 


7 


10.1 


.0810 


8 


12.5 


.0756 


9 


14.1 


.0699 


10 


15.6 


.0637 



TABLE III 

VALUES OF X AND A FOR VERTICAL 
n n 

MODES 1 THROUGH 6 



n 


A 

n 


X 

n 


1 


+.0237 


9. 18x10“5 


2 


-.0218 


4.66x10"^ 


3 


+.0193 


1.14x10“3 


4 


-.0169 


2.14x10“^ 


5 


+.0150 


3.47x10“3 


6 


-.0134 


5.15x10“^ 
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V. CONCLUSIONS 



The stated objective of this study was to model certain 
observed Arctic Ocean internal current patterns through an 
analytic study of the baroclinic response of the ocean to 
lead formation and refreezing. The problem was then divided 
into two major phases: (1) finding the baroclinic internal 

normal modes of the Arctic Ocean; and (2) investigation of 
the forcing of these modes by the effects of lead formation. 

The normal modes of response found as ei result of the 
theory developed in Section II showed a velocity pattern 
which was essentially stationary with no significant vertical 
velocity. However, horizontal velocities of the magnitude 
measured in the Arctic could exist within these stationary 
patterns. It was noted that the steepness of the pycnocline 
(measured by the parameter, c) had a direct and significant 
effect on the depth of maximum velocity, whereas the effect 
of a variation in 'layer depth was relatively minor. The 
apparent consistency of the cell structure seen in Figure 16, 
the theoretical vertical cross-section of horizontal veloci- 
ties, with the observed data lends strength to the plaus- 
ibility of the developed theory. 

The fact that a simple linear combination of the first 
two theoretical normal modes. Figure 14, qualitatively 
matches measured velocity profiles is of particular note. 

It shows the physics used to develop the theory to be 
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consistent with actual data, in that observed current 
patterns appear to be part of the set of allowed solutions. 

It has not yet been shown conclusively that the perturbations 
at a lead will force these inodes. 

The study done here is a first step which looked at the 
dynamics of a single Arctic lead. The .theory developed has 
been shown to be dynamically possible but the study has not 
yet looked at the energy involved in the system. It is 
possible that an interaction between leads or between lead(s) 
and the mean flow is needed to provide the energy required 
for the observed current magnitudes. 

The very small values of frequency, corresponding to 
periods of greater than 20 years, are attributable to the 
magnitude of 3 in the Arctic. The vorticity equation derived, 
eq. (38), is a non-linear time dependent parabolic equation. 
However, the results found in the initial value problem 
studied in Section IV of this thesis indicate that the time 
dependence of the equation is not significant since the solu- 
tions were essentially stationary. Under the same geometry 
the problem could be reformulated as an elliptic, linear 
equation forced with time dependent boundary conditions. 

An advantage of this type formulation v/ould be the ability 
to write a separable differential equation in cylindrical 
coordinates, which was not possible using the development of 
Section II. Motion approximating the vortex-like patterns 
observed by Coachman and Neviton [1973] could then be more 
easily described. 
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The plane wave solution for velocities developed in this 
thesis is a result of the formulation of the problem in 
rectangular coordinates with an infinite length lead. This 
worked v/ell in developing the theory and resulted in a real- 
istic vertical section of velocity structure. However, 
actual leads are of finite length which could be expected to 
induce curvature in the velocity patterns. 
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VI. 



RECOMMENDATIONS FOR FURTHER RFSEARCIl 



The consistency found with observed data warrants the 
continuation of this study using the theory developed in this 
thesis as a basis. The following recommendations for 
further research are made: 

1. the formulation of the initial value ’study as an 
elliptic problem, as discussed in section V; 

2. the use of vertical stretching caused by bathymetry, 
or forced convection at a lead, rather than the 3 effect as 
the driving mechanism in the developed vorticity equation; 

3. the investigation of the question of whether these 
velocity patterns have any influence on the stability of the 
larger scale mean flow, including the effects of baroclinic 
instabilities ; 

4. a study of the energetics of the system. 
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